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Abstract 

The  mechanical  properties  of  a  deformed  workpiece  are  sensitive  to  the  initial  mi¬ 
crostructure.  Often,  the  initial  microstructure  is  random  in  nature  and  location  specific. 
To  model  the  variability  of  properties  of  the  workpiece  induced  by  variability  in  the 
initial  microstructure,  one  needs  to  develop  a  reduced  order  stochastic  input  model  for 
the  initial  microstructure.  The  location-dependence  of  microstructures  dramatically  in¬ 
creases  the  dimensionality  of  the  stochastic  input  and  causes  the  “curse  of  dimension¬ 
ality”  in  a  stochastic  deformation  simulation.  To  quantify  and  capture  the  propagation 
of  uncertainty  in  multiscale  deformation  processes,  a  novel  data-driven  bi-orthogonal 
Karhunen-Loeve  (KL)  decomposition  strategy  is  introduced.  The  multiscale  random 
field  representing  random  microstructures  over  the  workpiece  is  decomposed  simul¬ 
taneously  into  a  few  modes  in  the  macroscale  and  mesoscale.  The  macro  modes  are 
further  expanded  through  a  second-level  KL  expansion  to  separate  the  random  and 
spatial  coordinates.  The  few  resulting  random  variables  are  mapped  to  the  uniform 
distribution  via  a  polynomial  chaos  (PC)  expansion.  As  a  result,  the  stochastic  input 
complexity  is  remarkably  simplified.  Sampling  from  the  reduced  random  space,  new 
microstructure  realizations  are  reconstructed.  By  collecting  the  properties  of  work- 
pieces  with  randomly  sampled  microstructures,  the  property  statistics  are  computed.  A 
high-dimensional  multiscale  disk  forging  example  of  FCC  nickel  is  presented  to  show 
the  merit  of  this  methodology,  and  the  effect  of  random  initial  crystallographic  texture 
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1.  Introduction 

For  polycrystalline  materials  (e.g.  metals,  alloys,  etc.),  macroscale  properties  are 
sensitive  to  the  underlying  microstructure.  In  this  work,  we  are  interested  in  model¬ 
ing  the  variability  of  properties  of  the  workpiece  in  a  deformation  process  induced  by 
variability  in  the  initial  microstructure.  To  accurately  evaluate  the  properties  of  a  given 
microstructure  at  a  given  material  point,  homogenization  strategy  over  a  statistical  vol¬ 
ume  element  (SVE)  [1,  2,  3]  has  been  widely  applied.  In  the  context  of  polycrystalline 
materials,  an  SVE  is  a  microstructure  containing  finite  number  of  grains  (see  Eig.  1) 
with  features  (grain  and  texture  distribution)  that  satisfy  certain  statistical  constraints. 
Given  a  number  of  correlated  realizations  of  the  microstructure,  a  stochastic  data-driven 
model  of  the  microstructure  is  produced  that  when  coupled  with  uncertainty  quantifi¬ 
cation  methods  (e.g.  Monte  Carlo  or  sparse  grid  methods)  can  be  used  to  compute 
the  distribution  of  properties  at  the  material  point.  However,  in  order  to  investigate  the 
variability  of  macroscopic  properties  in  a  workpiece  due  to  microstructure  randomness, 
we  need  to  exploit  the  stochastic  space  of  initial  microstructures  in  the  workpiece  and 
not  simply  at  a  material  point. 

Microstructure  uncertainty  at  a  material  point  has  been  extensively  studied  using  a 
variety  of  methods.  In  [4,  5]  the  principle  of  maximum  entropy  (MaxEnt)  was  used  to 
describe  the  microstructure  topology  of  binary  and  polycrystalline  materials.  A  set  of 
correlation  functions  or  grain  size  moments  were  given  as  the  prescribed  constraints. 
Realizations  of  microstructures  were  then  sampled  from  such  MaxEnt  distribution  and 
interrogated  using  appropriate  physical  model,  e.g.  a  crystal  plasticity  finite  element 
method  (CPEEM)  [6]  for  polycrystals.  The  Monte  Carlo  (MC)  method  was  adopted 
to  find  the  error-bars  of  effective  stress-strain  response  of  ECC  aluminum.  In  [7],  the 
orientation  distribution  function  (ODE)  was  adopted  to  describe  the  polycrystalline  mi- 
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crostmcture.  A  number  of  ODF  samples  were  given  as  the  input  data.  Karhunen-Loeve 
expansion  (KLE)  [8,  9]  was  utilized  to  reduce  the  input  complexity  and  facilitate  the 
high-dimensional  stochastic  simulation.  An  adaptive  version  of  sparse  grid  collocation 
strategy  [10,  11]  was  used  to  find  the  stress-strain  curve  with  error  bars  and  the  convex 
hull  of  elastic  modulus  of  FCC  aluminum  after  deformation.  Mechanical  response  vari¬ 
ability  and  thermal  properties  due  to  both  orientation  and  grain  size  uncertainties  were 
studied  in  [12,  13].  A  nonlinear  model  reduction  technique  based  on  manifold  learn¬ 
ing  [14]  has  been  introduced  to  find  the  surrogate  space  of  the  grain  size  feature  while 
grain  orientations  were  reduced  by  KLE.  Critical  stress  distribution  after  deformation 
was  constructed  for  FCC  nickel  [12]  and  effective  thermal  conductivity  distribution 
was  explored  for  crystals  with  ortho-symmetry  [13].  Recently,  variability  of  fatigue 
resistance,  measured  by  strain-based  fatigue  indicator  parameters  (FIPs)  [15]  of  two- 
phase  nickel-based  superalloys  was  studied  with  the  assistance  of  principal  component 
analysis  (PC A)  [16].  Distributions  of  FIPs,  as  well  as  their  convex  hulls  showing  the 
extreme  values,  of  microstructures  sharing  identical  statistical  features  with  given  sam¬ 
ples  under  cyclic  loading  were  extracted.  Convergence  with  increasing  dimensionality 
of  the  reduced-order  representation  was  also  shown. 

To  quantify  the  effect  of  microstructures  on  macro-properties  and  probe  the  un¬ 
certainty  propagation  through  different  length  scales,  a  multiscale  simulator  needs  to 
be  adopted.  Each  point  of  the  workpiece  is  associated  with  a  microstructure  in  the 
mesoscale,  the  deformation  of  which  is  controlled  by  the  local  deformation  gradient 
estimated  in  the  macroscale.  Mechanical  properties/response  of  the  point  are  evalu¬ 
ated  via  proper  (e.g.  crystal  plasticity)  constitutive  model  applied  on  the  deformed 
microstructure.  Since  microstructures  are  random,  properties  of  the  workpiece  are  also 
random.  In  general,  microstructures  are  location-specific  (meaning  that  microstruc¬ 
tures  associated  with  different  spatial  points  may  have  different  distributions)  [17].  As 
a  result,  the  stochastic  input  to  the  probabilistic  multiscale  forging  simulation  will  be 
extremely  high  dimensional,  which  prevents  one  from  quantifying  the  uncertainties  of 
interested  properties.  This  problem  is  usually  referred  to  as  the  “curse  of  dimensional¬ 
ity”.  Conventional  model  reduction  schemes  that  only  locally  decompose  input  com¬ 
plexity  within  a  single  scale  and  cannot  see  the  correlation  between  macroscale  points 
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are  not  sufficient  for  reducing  the  multiscale  stochastic  input.  To  this  end,  we  intro¬ 
duced  a  bi-orthogonal  KLE  strategy  [18, 19],  which  decomposes  the  multiscale  random 
field  into  a  few  modes  in  both  the  macro-  and  meso-scales  [20].  Mean  and  standard 
deviation  of  elastic  moduli,  i.e.  Young’s  modulus,  shear  modulus,  and  bulk  modulus, 
over  the  product  were  investigated  given  texture  (orientation  distribution)  uncertainty 
of  FCC  copper  microstructures.  However,  this  earlier  work  limited  its  stochastic  in¬ 
put  to  two  prescribed  random  variables.  In  this  work,  we  will  build  the  bi-orthogonal 
KLE  on  the  basis  of  a  given  set  of  microstructure  data.  Different  types  of  macroscale 
inner  products  are  discussed.  A  second-level  KLE  is  conducted  to  further  reduce  the 
dimensionality  of  the  stochastic  space  after  bi-orthogonal  decomposition.  The  optimal 
dimensionality  of  the  final  reduced-order  space  will  be  determined  based  on  the  energy 
proportion  captured  by  the  principal  components  in  the  two-step  decomposition.  A 
non-intrusive  strategy  is  used  to  project  the  reduced  random  variables  to  the  space  of 
random  variables  with  known  probability  distributions.  Low  order  statistics  of  equiv¬ 
alent  stress,  strain,  and  strength  fields  of  disks  after  forging  are  studied  by  repeatedly 
calling  the  deterministic  solver  using  microstructures  sampled  in  the  reduced  space. 
We  use  Monte  Carlo  (MC)  sampling  to  construct  the  stochastic  solution. 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2,  the  representation  of 
microstructure  and  the  multiscale  forging  simulator  are  introduced.  Theory  and  formu¬ 
lation  of  the  bi-orthogonal  KL  decomposition,  followed  by  a  second-level  KLE  to  fur¬ 
ther  decompose  the  resulting  spatial-random  coupled  modes,  will  be  enunciated  in  Sec¬ 
tion  3.  The  reduced  microstructure  representations  are  mapped  to  a  multi-dimensional 
uniform  distribution  via  polynomial  chaos  expansion  (PCE)  so  that  sampling  of  mi¬ 
crostructures  becomes  efficient.  The  detail  of  this  procedure  is  presented  in  Section  4. 
Problem  definition  and  generation  of  initial  samples  are  described  in  Section  5.  Sec¬ 
tion  6  shows  an  example  of  stochastic  multiscale  modeling  of  nickel  disk  forging.  Prop¬ 
erties/response  variability  of  the  product  will  be  investigated.  Finally,  conclusions  and 
discussion  are  given  in  Section  7. 
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2.  Multiscale  Modeling  of  a  Forging  Process 


To  study  the  variability  of  mechanical  properties  induced  by  microstructure  un¬ 
certainty  in  forging  disks,  a  multiscale  framework  which  couples  finite  element  (FE) 
large  deformation  simulator  with  crystal  plasticity  constitutive  model  is  introduced  as 
the  deterministic  solver.  Each  point  in  the  macroscale  is  linked  to  a  mesoscale  poly¬ 
crystalline  microstructure  described  by  its  grain  size  and  orientation  features.  The  me¬ 
chanical  properties  and  response  of  the  material  under  deformation  are  evaluated  in  the 
microstructure  domain  and  returned  to  the  workpiece. 

2.7.  Multiscale  Forging  Solver 

An  updated  Lagrangian  implicit  EE  model  for  the  analysis  of  large  deformation 
forging  processes  is  employed  for  the  multiscale  simulation.  This  model,  seamlessly 
integrating  kinematic,  contact,  and  constitutive  modules,  is  suitable  for  forging  problem 
of  any  material  subjected  to  various  die  shapes.  A  crystal  plasticity  model  is  used  in  this 
work.  Each  point  (effectively.  Gauss  integration  point)  of  the  macroscale  workpiece  is 
represented  by  a  polycrystalline  microstructure  whose  deformation  is  controlled  by  the 
local  deformation  gradient  F,  while  the  mechanical  properties/response  (e.g.  effective 
stress,  strain,  strength,  etc.)  at  that  point  are  evaluated  in  the  mesoscale  and  returned 
to  the  workpiece  after  homogenization  [21].  The  multiscale  framework  is  summarized 
below. 

Macroscale:  Let  X  be  a  material  particle  in  the  initial  configuration  So  and  x  = 
x(X,  f„+i)  be  its  location  at  time  tn+\  -  The  total  deformation  gradient  defined  as 

F(X,  f„+i)  =  Vox(X,  Wi)  =  (1) 

can  be  expressed  as  the  product  of  F„  at  the  previous  time  step  t  —  tn  and  the  relative 
deformation  gradient  F^; 

F  =  F,F„.  (2) 

The  equilibrium  equation  at  f  =  f„+i  is  expressed  in  the  reference  configuration  S„ 
as 

V„  ■  <P,)a  +  f.  =  0,  (3) 
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where  V„  denotes  the  divergence  in  S„.  can  be  represented  as 


fr  =  det  Frb, 


(4) 


where  b  is  the  body  force  defined  on  the  current  configuration.  The  material  behavior 
is  obtained  from  the  deformation  of  the  microstructure  through  homogenization.  The 
function  (■)/,  denotes  the  homogenized  property  over  the  microstructure.  Therefore, 
{Pr)h  is  the  homogenized  first  Piola-Kirchhoff  (PK)  stress  expressed  per  unit  area  of 
Sn  and  given  as: 

(P,),  =  (detF^TFr^),  =  detF,<T)ftF;^  (5) 

where  F^  is  the  mesoscale  relative  deformation  gradient,  which  equals  the  macroscopic 
Fr  under  Taylor  hypothesis  for  macro-meso  linking.  (T)/,  is  the  homogenized  Cauchy 
stress  defined  as  the  volume  average  of  the  mesoscale  stress  T: 

<T),  =T  =  i  r  TdV.  (6) 

An  incremental  quasi-static  problem  is  solved  to  determine  the  displacement  field 
u(x„,  tn+i)  that  satisfies  Eq.  (3)  and  boundary  conditions.  The  solution  of  the  deforma¬ 
tion  problem  proceeds  incrementally  in  time  starting  from  the  initial  configuration  Sq. 
The  weak  form  of  the  governing  equation  in  the  presence  of  die  contact  can  be  written 
as: 

G(u„+i,u)  =  Gi(u„+i,u)  -hGc(u„+i,u)  =  0,  (7) 


where  the  first  term  is  the  virtual  work  of  the  workpiece: 

G/,(u„+i,u)  =  r  (P}h--^dV 

-  (  f  t  ■  QdS  +  f  b  ■  udV  I , 

Js„+i 


(8) 


and  the  second  term  is  the  contact  virtual  work: 


Gc(u„+1,  u)  =  I  (-w  ■  U  -H  tcT-  ■  u)ds,  (9) 

Jas‘„ 

where  u  is  the  virtual  displacement;  t  is  the  traction  and  b  is  the  body  force  of  the  bulk; 
c  dS„  is  the  surface  corresponding  to  regions  of  the  body  that  may  potentially 
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come  into  contact  with  the  die.  tcA^  and  tej-  are  the  normal  and  tangential  tractions  at 
the  die  due  to  contact  (and  friction).  The  contact  work  is  calculated  in  the  reference 
configuration  S„.  Newton-Raphson  iterative  scheme  along  with  a  line  search  method  is 
used  to  solve  this  non-linear  system.  The  linearization  of  the  weak  form  at  the  (k-H  l)-th 
iteration  becomes 


G(uri„u)  + 


dG 


du 


(k) 

n+1 


■(U 


(t:+l) 


-  U®  ^  - 
n+l  ^n+1 


)  =  o. 


(10) 


The  increment  of  internal  work  Gi,  is  computed  by 


dGh^  f  d(P,h-^dV.  (11) 

Js„  dx„ 

The  linearization  of  the  homogenized  PK-I  stress  is 


d(Pr}k  =  detF,(fr(ufF,F;')<T)ft 

-  (Th(dF,F;Y  +  (dTh)F;^.  (12) 


The  homogenized  Cauchy  stress  (T)/,  and  its  gradient  with  respect  to  dF^  are  evaluated 
using  the  constitutive  model  in  the  mesoscale. 

Mesoscale:  As  mentioned  earlier,  each  Gauss  point  of  the  workpiece  in  the  macroscale 
corresponds  to  a  polycrystalline  microstructure,  which  is  described  by  the  sizes  and  ori¬ 
entations  of  its  constituent  grains.  According  to  the  Taylor  hypothesis,  all  grains  of  the 
same  microstructure  are  subject  to  the  same  deformation  gradient  F,  which  is  identical 
to  the  local  deformation  gradient  at  that  point  of  the  workpiece.  The  mechanical  re¬ 
sponse  of  each  grain  is  computed  using  the  crystal  plasticity  constitutive  model  [6,  22] 
and  averaged  over  the  microstructure  to  represent  the  corresponding  quantity  in  the 
macroscale.  In  this  paper,  we  are  interested  in  FCC  nickel  polycrystals.  For  details  of 
the  constitutive  model,  the  reader  is  referred  to  [23].  We  here  only  briefly  review  the 
main  algorithm. 

The  total  deformation  gradient  F  is  multiplicatively  decomposed  into  elastic  and 
plastic  parts 

F  =  F'F^,  (13) 

where  F^  is  the  elastic  deformation  gradient  and  F^  is  the  plastic  counterpart  with 
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detF^  being  1.  The  PK-II  stress  is  computed  by 

t  =  -  I),  (14) 

where  Xf  is  the  fourth-order  elasticity  tensor  represented  in  the  sample  coordinate  sys¬ 
tem  and  I  is  the  second-order  identity  tensor. 

The  resolved  shear  stress  on  slip  system  a  is  calculated  as 

=  t  ;  (15) 

where  s  ®  is  the  Schmid  tensor  defined  by  the  tensor  product  of  the 
slip  direction  and  slip  normal  of  slip  system  a.  If  is  larger  than  the  slip 
resistance  plastic  flow  on  the  a-th  slip  system  occurs.  The  incremental  plastic 
shear  strain  is  therefore  obtained  according  to  flow  rule  and  utilized  to  update  the 
slip  resistance,  deformation  gradient,  etc.  The  Taylor  strain  hardening  law  based  on 
dislocation  density  p  is  adopted  for  the  homogeneous  evolution  of  slip  resistance  k. 

K  —  kq  +  atpb  yfp,  (16) 


where 

or  '  ^  / 

and  at  is  a  constant  representing  an  average  of  the  junction  strength  over  all  existing 
dislocation  configurations,  p  is  the  shear  strength,  and  b  is  the  magnitude  of  the  Burgers 
vector.  The  first  term  in  Eq.  (17)  represents  a  geometric  storage  due  to  lattice  incompat¬ 
ibility,  describing  the  grain  boundary  hardening.  Lg  is  the  grain  size  parameter  [  1 2,  24] . 
The  second  term  describes  storage  through  a  statistical  measure  of  forest  dislocation, 
describing  the  dislocation  interaction  hardening  inside  grains.  The  last  term  represents 
a  dynamic  recovery  rate  that  renders  dislocation  segments  inactive  as  they  rearrange 
themselves  [22].  The  parameters  ki  and  k2  are  given  as 


h 


200  ^  _  20q 

Utpb  ’  Ks-  Kq’ 


(18) 


where  0q  is  the  initial  hardening  rate,  kq  and  Ks  are  the  initial  yield  stress  and  saturation 
strength,  respectively.  Substituting  Eqs.  (17)  and  (18)  into  Eq.  (16),  we  can  derive  at 


the  hardening  rate  as 


2L„{k  -  kq) 


+  Sq 


Ks-  K 
Ks  -  Kq 


Therefore,  the  slip  resistance  at  f  =  f„+i  is  updated  as 

Kn+\  —  Kn+  klS.  t. 


(19) 


(20) 


An  iterative  algorithm  is  designed  to  solve  this  system  of  nonlinear  equations.  Macro¬ 
scopic  quantities,  such  as  stress  and  strain,  are  computed  as  the  volume-average  of  the 
mesoscale  values  over  all  grains  (e.g.  Eq.  (6)).  The  macroscopic  von-Mises  equivalent 
stress  and  equivalent  strain  are  calculated  in  the  form  of 


^eff-  x  It'-T', 


where  T'  is  the  deviatoric  part  of  the  homogenized  Cauchy  stress  T,  and 


^eff 


■Iff 


Ddt, 


(21) 


(22) 


in  which  D  is  the  average  rate  of  deformation.  Grain  orientations  evolve  accordingly 
due  to  elastic  deformation; 


m 


(a)  _ 


F"m 


(a) 


n(“)  = 


(23) 


The  equivalent  strength  is  evaluated  as  the  average  slip  resistance  of  all  slip  systems 
of  all  grains  in  the  microstructure: 


f^slip 


Xa) 


(24) 


The  homogenized  properties  and  response  are  returned  to  the  macroscale  for  up¬ 
dating  the  deformation  and  response  fields  of  the  workpiece.  The  verification  of  this 
particular  crystal  plasticity  constitutive  model  is  discussed  in  [12]. 


2.2.  Microstructure  Representation 

In  the  current  work,  the  Taylor  hypothesis  is  adopted  for  the  macro-mesoscale 
linking.  Accordingly,  a  microstructure  (Fig.  1(a))  can  be  simply  described  by  an  ar¬ 
ray  of  topological  and  orientational  features  of  constituent  grains  (Fig.  1(b)).  For  a 
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microstmcture  (e.g.  FCC  nickel)  composed  of  M  grains,  the  first  M  components  of 
the  feature  array  are  sizes  of  individual  grains  sorted  in  ascending  order  and  the  rest 
3M  components  are  the  corresponding  orientations  described  by  Rodrigues  parame¬ 
ters  [25],  an  axis-angle  representation  that  consists  of  three  components  dehned  in 
Eq.  (25): 


r  =  w  tan  - , 
2 


(25) 


where  r  =  {r\,r2,  r3)  are  the  three  Rodrigues  components;  w  =  {wi,  W2,  W3)  gives  the 
direction  cosines  of  the  rotation  axis  with  respect  to  microstructure  coordinates;  and  6 
is  the  rotation  angle. 


Figure  1 :  (a)  A  3D  polycrystalline  microstructure  with  20  grains,  (b)  The  descriptor  of  the  microstructure. 
The  first  20  components  are  the  sizes  of  grains,  and  the  last  60  components  are  the  Rodrigues  parameters 
representing  grain  orientations. 


The  poly  cry  stalline  microstructure  is  usually  a  high-dimensional  random  held,  which 
makes  the  stochastic  simulation  intractable.  For  example,  a  20-dimensional  vector  is 
needed  to  store  the  grain  size  feature  of  a  microstructure  containing  20  grains.  The 
vector  will  end  up  to  be  80-dimensional  when  Rodrigues  parametrization  orientations 
are  considered  as  well.  The  situation  is  even  worse  when  the  random  held  is  also  a 
function  of  spatial  location  x.  If  the  correlation  between  microstructures  at  different 
points  on  macroscale  is  not  explored,  the  dimensionality  of  the  random  held  explodes. 
For  a  random  2D  workpiece  discretized  by  Uei  quadrilateral  elements  each  of  which 
has  riint  Gauss  points,  the  total  dimensionality  of  the  microstructure  descriptor  ends  up 
to  be  AMneitiint,  where  M  is  the  number  of  grains  in  the  microstructure.  In  the  current 
work,  we  assume  the  random  source  is  only  the  grain  orientations  while  the  grain  sizes 
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are  fixed.  The  total  dimensionality  of  the  microstructure  stochastic  space  is  therefore 
hMueiriint,  which  is  still  very  large.  This  is  referred  to  as  the  “curse  of  dimensional¬ 
ity”.  A  reduced-order  surrogate  microstructure  model  of  the  location  dependent  ran¬ 
dom  microstructure  is  needed.  By  sampling  from  the  low-dimensional  surrogate  space, 
uncertainty  quantification  of  the  product  properties  driven  by  random  microstructures 
becomes  computationally  feasible. 

3.  Two-Step  Karhunen-Loeve  Decomposition  of  the  Multiscale  Random  Microstruc¬ 
ture  Field 

In  our  previous  works  on  quantifying  uncertainties  in  materials,  a  set  of  linear 
and  non-linear  model  reduction  techniques  were  developed  to  facilitate  the  solution 
of  stochastic  partial  differential  equations  (SPDEs)  describing  physical  processes  in 
random  media  [14,  26,  27,  28].  These  methods  have  been  successfully  applied  to  re¬ 
duce  the  dimensionality  of  random  microstructures  at  a  given  material  point.  However, 
these  techniques  cannot  be  applied  to  location-dependent  microstructures  (Fig.  2),  e.g. 
reducing  a  random  microstructure  field  of  the  form  A(x,  s,  (o),  where  s  is  a  mesoscale 
coordinate  and  x  the  macroscopic  coordinates.  A  bi-orthogonal  decomposition  based 
on  KL  expansion  was  introduced  to  address  this  problem  in  [20].  This  algorithm  was 
originally  developed  for  temporal-spatial  coupled  problems  [18],  in  which  the  random 
field  ,  A(x,  f,  cj),  is  a  function  of  both  time  t  and  spatial  coordinates  x.  Following  cer¬ 
tain  modification,  we  have  applied  this  idea  to  multiscale  forging  problem  in  [20].  In 
the  current  work,  we  rebuild  the  bi-orthogonal  KFF  strictly  on  the  basis  of  given  data 
and  generalize  the  model  to  the  scenario  that  the  information  of  the  inherent  controlling 
random  variables  is  not  known  beforehand.  The  stochastic  random  field  A(x,  s,  (o)  is 
decomposed  to  a  set  of  mesoscale  modes  'I'(s)  and  macro  (or  spatial)-random  coupled 
modes  0(x,  6>).  A  second-level  KFF  is  introduced  to  further  separate  random  variables 
(o  from  the  macroscale  (spatial)  modes.  The  optimal  dimensionality  of  the  reduced- 
order  space  will  be  determined  based  on  the  number  of  principal  components  preserved 
in  the  two-step  KFF.  Different  selections  of  the  inner  product  in  the  macroscale  are  dis¬ 
cussed  to  determine  which  one  is  more  appropriate  for  this  problem.  The  algorithm  of 
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bi-orthogonal  decomposition  (BOD)  is  summarized  below. 


Figure  2:  Microstracture  dependence  on  spatial  location.  At  different  locations,  the  microstructures  may 
have  different  features  due  to  preprocessing. 


3.1.  Bi-Orthogonal  Karhunen-Loeve  Decomposition 

Assume  a  random  field  A  defined  on  a  probability  space  (Q,  F,  p) 

A(x,  s,  to) :  .S  X  Ad  X  Q  ^  K,  (26) 

where  S  is  the  macroscale  spatial  domain,  M  is  the  microstructure  space,  Q.  is  the  set 
of  elementary  events  and  tu  e  Q  is  a  random  field  that  determines  the  uncertainty  of  A 
(we  call  it  controlling  random  variable).  A  =  {A',  /  =  1, . . . ,  m)  can  be  the  combination 
of  m  independent  features,  e.g.,  grain  size,  texture,  etc.,  for  describing  a  microstructure. 
In  the  example  shown  later,  m  is  chosen  to  be  1  as  crystallographic  texture  is  the  only 
random  source.  The  formulation  presented  in  this  section  takes  the  general  situation 
where  A  is  the  combination  of  various  features.  One  can  use  the  idea  of  Karhunen- 
Loeve  expansion  to  project  this  field  to  a  set  of  bi-orthogonal  bases  in  the  form  of 

A(x,  s,  (o)  -  A(x,  s)  -H  A(x,  s,  OJ) 

=  A(x,  s)  +  J  ^T'f)(s)Of >(x,  oj),  ill) 

i=\ 

in  which  p®  are  eigenvalues  of  the  eigenvalue  problem  to  be  defined  later  on,  are 
mesoscale  modes  strongly  orthogonal  in  the  meso-space,  and  O,-  =  {(!)[**’*, . . .  (p®’'") 
are  spatial  modes  weakly  orthogonal  in  macro-space  for  all  the  m  features.  The  su¬ 
perscript  (h)  =  1,2,3  denotes  that  we  can  construct  the  above  expansion  according  to 
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certain  type  of  inner  products  in  the  macroscale.  A  is  the  mean  field  defined  by 

A  =  (A)  I  A(x,  s,  (28) 

Jn 

where  p(io)  is  the  multivariate  joint  probability  density  of  to.  The  definition  of  A  is  not 
unique.  For  example,  in  [19],  this  mean  field  is  also  averaged  over  the  temporal  domain 
(analogous  to  the  mesoscale  space  in  the  current  setup),  leaving  A  a  function  of  spatial 
coordinates:  A  =  A(x).  We  here  do  not  average  it  over  the  meso  domain  to  be  con¬ 
sistent  with  the  previous  work  [20].  In  practice,  the  random  field  A(x,  s,  a>)  is  usually 
given  by  N  realizations  (A;(x,  s,  (y,)]^j.  As  a  consequence,  the  mean  field  is  computed 
as  the  average  of  all  given  samples.  For  model  reduction,  the  sum  in  Eq.  (27)  is  usually 
approximated  by  the  first  finite  number  of,  say  d,  principal  components  (modes)  that 
capture  most  of  the  energy. 

We  denote  by  (, )  the  inner  product  in  the  microstructure  domain  and  hy  {,}fi  (h  - 
1, 2,  or  3)  different  types  of  inner  products  in  the  spatial  domain  [19],  respectively,  and 
obtain 

r  'I'‘'‘\s)'T®(s)t/s,  (29) 

^  Jm  ^ 

and 
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(h) 

The  mesoscale  modes  are  strongly  orthogonal  to  each  other  and  the  macroscale 
modes  O®  are  weakly,  in  an  average  sense,  orthogonal. 

By  minimizing  the  distance  (based  on  the  norm  defined  in  Eq.  (30))  between  the 
Karhunen-Loeve  expansion  and  the  original  random  field,  one  ends  up  with 


{A,of>U, 


and  from  the  orthogonality  condition,  we  obtain 


w)  = 


1 


f 

Jm 


Jh)  Jm 


A(x,  s,  (oyvf^{s)ds. 


(33) 


(34) 


These  last  two  Eqs.  (33)  and  (34)  lead  to  the  following  eigenvalue  problem 

C®(s,s)'I'f\s)t/s, 


pW'Pf>(s)-  r 

J  A 


M 

where  the  covariance  matrix  is 

s)  =  {A(x,  s,  (jj),  A(x,  s,  (Ji>)]h- 
In  discrete  form,  the  covariance  can  be  written  as 


^el  f^int 


CW(s,s)  - 

iii  =  l  im  =  ^  V 


N 


\ 


1  ^  \ 
N  A-i  ‘m’  ■' 


1 


N  riel  mm 


c'^^(s,s)  =  -^^^a/x;2,s)aJ(x;2,s)b^j4I, 

j=l  i^=l  i^  =  l 

C<2)(s,s)  =  C<iVs,s)-C<°Vs,s), 


(35) 


(36) 


(37) 


for  different  definition  of  macroscale  inner  products.  In  these  equations,  N  is  the  num¬ 
ber  of  realizations,  Wg;  is  the  number  of  elements  in  macroscale  (when  a  EEM  dis¬ 
cretization  of  the  domain  is  created),  n,„,  is  the  number  of  integration  points  in  each  el¬ 
ement,  \JiJ  is  the  Jacobian  determinant  of  the  element  is  the  integration  weight 

associated  with  the  integration  point  and  A  is  a  matrix  containing  microstructural 
features  corresponding  to  integration  points,  and  x'"  represents  global  coordinates  of 
the  integration  point  of  element  /„  in  the  macroscale. 
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The  three  types  of  covariance  matrix  defined  in  Eq.  (37)  (or  equivalently,  the  in¬ 
ner  products  in  the  macroscale  Eq.  (30))  lead  to  constructions  of  different  expansions 
Eq.  (27)  according  to  different  optimization  of  the  random  field  A:  through  /z  =  0,  we 
minimize  the  mean  error;  through  h  -  1 ,  we  minimize  the  second-order  moment  of 
the  error  (the  Euclidian  distance  between  the  KL  expansion  and  the  random  field);  and 
through  h  —  2,  we  minimize  the  standard  deviation  error.  The  discussion  of  the  three 
types  will  be  further  enunciated  later  with  a  practical  example. 

3.2.  A  Second-Level  KLE 

With  the  bi-orthogonal  KL  expansion,  we  decompose  the  random  field  A(x,  s,  ai) 
into  a  set  of  mesoscale  modes  {T',)  and  spatial-random  coupled  modes  (O,  ).  The  dimen¬ 
sionality  of  the  original  stochastic  microstructure  space  can  be  reduced  by  truncating 
the  bi-orthogonal  KLE  (Eq.  (27))  to  d  terms,  which  capture  most  of  the  “energy”  of  the 
given  samples.  The  truncated  bi-orthogonal  KLE  is  given  by  Eq.  (38). 

A(x,  s,  (O)  A(x,  s)  H-  ^  ^/^T'f\s)OfVx,  (O).  (38) 

1=1 

The  mesoscale  modes  are  known  as  the  eigenvectors  of  the  covariance  matrix  and 
only  depend  on  the  mesoscale  coordinates  s.  The  macroscale  modes,  on  the  other 
hand,  couple  the  random  source  (o  with  spatial  coordinates  x,  and  they  are  still  high¬ 
dimensional  random  fields.  To  further  reduce  the  dimensionality,  the  spatial-random 
coupled  modes  resulting  from  the  bi-orthogonal  KL  expansion  need  further  decompo¬ 
sition.  A  conventional  strategy  is  to  separate  the  random  variables  a>  from  the  spatial 
coordinates  x  using  a  polynomial  chaos  expansion  (Eq.  (39)), 

0f>(x,  tu)  =  ^  7W(x)Tf(^(6>)),  (39) 

j 

where  T®(^((y))  are  orthogonal  polynomials  of  random  variables  ^(a>),  which  usually 
follow  well-known  probability  distributions.  In  the  current  problem,  however,  the  com- 
putation  of  location-specific  coefficients  yjj  ("s.)  is  very  complicated  as  O;  (x,  cl>)  are 
random  fields  in  terms  of  spatial  coordinates  x  and  are  known  only  in  the  form  of  N 
(equal  to  the  number  of  initial  samples)  finite  number  of  realizations  derived  from  the 
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bi-orthogonal  KLE.  The  PCE  cannot  be  applied  in  a  straightforward  manner  to  these 
macro-random  coupled  modes. 

We  hereby,  propose  two  reasonable  assumptions  to  simplify  the  problem  and  em¬ 
ploy  a  second-level  KLE  to  decompose  the  random  variables  from  spatial  coordinates. 
The  assumptions  are: 

•  The  inherent  controlling  random  variables  w  can  be  separated  from  the  mesoscale 
and  macroscale  coordinates  (s,  x),  as  well  as  the  features  of  the  microstructure 
(i.e.  the  randomness  is  independent  of  the  microstructure  features,  e.g.  texture 
and  grain  size  have  the  same  random  source). 

(h) 

•  The  macroscale  modes,  are  independent  from  each  other. 

The  first  assumption  is  natural  since  it  is  the  fundamental  of  the  decomposition. 
The  second  assumption  is  a  strong  assumption  for  arbitrary  stochastic  processes,  since 
we  only  have  the  weak  orthogonality  condition  between  macroscale  modes  (Eq.  (32)). 
However,  it  is  important  for  the  further  decomposition  of  the  spatial  modes  0,(x,  ai) 
and  we  will  see  later  on  in  this  paper  that  this  assumption  leads  to  accurate  results. 

Having  the  above  two  assumptions,  we  can  next  perform  a  second-level  KLE  on 
each  macro-random  coupled  mode  O,.  Omitting  the  subscript  (h),  the  /-th  macro-mode 
O  ,(x,  a>i)  can  be  expanded  by 

r  ——  N 

0,(x,  iOi)  ^  0,(x)  +  £  7^'^/(x)^/(tu,),  ^  (40) 

j=l  i=l 

where  and  ij/^.  are  the  j-th  eigenvalue  and  eigenvector,  respectively,  of  the  covariance 
matrix, 

=  -*'■)’  (41) 

^=1 

and  r,  is  the  number  of  largest  eigenvalues  i'ijVjLi  that  capture  most  of  the  energy  of 
the  samples.  where  i  -  1, ...  ,c/,  are  the  reduced  representations  of  the  original 

multiscale  microstructures. 

In  this  way,  the  macro-random  coupled  modes  {<!>,  )  are  decomposed  into  basis  func¬ 
tions  depending  on  the  macroscale  coordinates  x  and  uncorrelated  random  variables 
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{<!>{]■  Note  that  the  dimensionality  of  the  stochastic  space  is  reduced  for  the  second 
time.  The  dimensionality  of  the  final  reduced  random  space  of  microstructures  over 
the  macroscale  is  the  sum  of  the  principal  dimensions  that  are  preserved  for  represent¬ 
ing  macro-random  modes: 

d 

(42) 

i=l 

where  d  is  the  truncated  dimensionality  of  the  bi-orthogonal  KLE  (Eq.  (38)).  The 
reduced  stochastic  space  can  then  be  constructed  and  equivalently  mapped  to  well- 
known  probability  distributions  through  polynomial  chaos  expansion  introduced  in  the 
next  section. 

4.  Polynomial  Chaos  Expansion  of  Stochastic  Reduced-Order  Model 

After  obtaining  reduced  representations  ({^/)^‘_j,  where  i  -  1, . . . ,  t/)  of  the  multi¬ 
scale  microstructure  samples,  we  next  need  to  construct  the  reduced-order  space  and 
map  it  to  a  probabilistic  space  from  which  new  samples  are  easy  to  draw.  Any  sample 
generated  in  the  low-dimensional  surrogate  space  can  be  recovered  to  spatial  modes  O 
through  Eq.  (40)  and  further  reconstructed  to  a  new  multiscale  microstructure  realiza¬ 
tion  in  the  original  input  space.  Polynomial  chaos  expansion  [9,  29,  30]  is  therefore 
introduced  to  represent  as  a  function  of  Gaussian  or  uniform  random  variables  As 
mentioned  before,  components  of  are  uncorrelated  but  not  necessarily  indepen¬ 

dent.  Although  Rosenblatt  transformation  [31]  can  be  used  to  decompose  the  problem 
to  a  set  of  independent  random  variables,  this  is  computationally  expensive,  especially 
for  high-dimensional  problems.  In  this  work,  we  assume  the  independence  between 
the  components  of  It  has  been  shown  in  various  applications  [32,  33,  28,  16] 

that  this  assumption  gives  rather  accurate  results. 

Eollowing  the  independence  assumption  of  <p^.,  each  of  them  can  be  expanded  on  to 
an  one-dimensional  polynomial  chaos  (PC)  basis  of  degree  p: 

p 

k=0 

where  the  are  i.i.d.  random  variables.  The  random  basis  functions  are  chosen 

according  to  the  type  of  random  variable  that  has  been  used  to  describe  the  random 
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input.  For  example,  if  Gaussian  random  variables  are  chosen  then  the  Askey  based 
orthogonal  polynomials  {T^)  are  chosen  to  be  Hermite  polynomials;  if  are  chosen  to 
be  uniform  random  variables,  then  {T*)  must  be  Legendre  polynomials  [29]. 

Uniform-Legendre  format  is  taken  for  the  projection  of  the  reduced-order  random 
variables  as  it  lends  very  close  reconstruction  of  the  PDFs  of  reduced  variables  (see 
Section  6).  The  PC  coefficients  are  computed  by 


(44) 


When  Uniform-Legendre  is  chosen,  Eq.  (44)  becomes 

rf  ^  ^  J  k^o,...,p.  (45) 

A  proper  method  is  needed  to  evaluate  these  integrals.  It  is  noted  that  the  random 
variable  does  not  belong  to  the  same  stochastic  space  as  and  we  only  have  N 
realizations  of  (fij.  The  distribution  of  is  not  known.  A  non-linear  mapping  F  :  ^ 

is  thus  needed  which  preserves  properties  such  that  and  have  the  same 
distribution.  A  non-intrusive  projection  based  on  empirical  cumulative  distribution 
functions  (CDFs)  of  samples  developed  in  [33]  is  utilized  to  build  the  map.  The  integral 
in  Eq.  (45)  is  then  computed  using  Gauss  quadrature. 

The  non-linear  mapping  F  ;  ^  ^  can  be  defined  as  shown  below  for  each  : 


=  F/(0, 


(46) 


where  F^j  and  F^j  denote  the  CDFs  of  and  respectively.  Here,  the  equalities, 
is  interpreted  in  the  sense  of  distribution  such  that  the  probability  density  functions 
(PDFs)  of  random  variables  on  both  sides  are  equal.  The  marginal  CDF  of  samples 
can  be  evaluated  numerically  from  the  available  data.  Kernel  density  estimation  is  used 
to  construct  the  empirical  CDF  of  (f>j.  Let  be  the  N  samples  of  obtained 

from  Eq.  (40).  The  marginal  PDE  of  is  then; 


(47) 
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The  marginal  CDF  of  is  obtained  by  integrating  Eq.  (47)  and  the  inverse  CDF  can  be 
computed  as  well.  Having  the  map  r(,  the  coefficients  yf  are  subsequently  computed 
via  Gauss  quadrature. 

After  mapping  the  reduced  space  to  uniform  distribution  random  variables,  the 
Monte  Carlo  method  can  be  used  to  sample  new  microstructure  realizations. 


5.  Problem  Definition  and  Initial  Data  Set 


The  problem  of  interest  in  this  work  is  the  variability  of  mechanical  properties 
of  forged  disks  due  to  (initial)  microstructure  uncertainties.  The  very  first  task  is  to 
generate  initial  correlated  random  microstructures  of  workpieces  to  be  used  as  the  input 
to  the  stochastic  simulation.  The  initial  workpiece  samples  are  a  set  of  cylindrical 
(rectangular  in  a  2D  axisymmetric  representation)  ingots,  each  point  of  which  is  finked 
to  a  distinct  microstructure  (Fig.  3). 


Figure  3:  Input  to  the  multiscale  deformation  simulator. 


We  generate  realizations  of  tfie  preforms  witfi  corTelated  microstructure  using  a  pre¬ 
processing  tfiat  deforms  a  set  of  raw  ingots  witfi  random  surfaces  into  regular  cylinders. 
The  simulations  are  done  within  a  2D  axisymmetric  Lagrangian  FE  framework.  The 
ingots  are  discretized  by  10x6  quadrilateral  elements,  each  of  which  contains  4  Gauss 
points  for  the  integration  in  the  element  domain.  Each  Gauss  point  is  linked  to  a  mi¬ 
crostructure  consisting  of  20  grains.  One  thousand  raw  ingots  whose  upper  surfaces 
(Fig.  4)  are  represented  by  degree  6  Bezier  curves  (Eq.  (48))  are  firstly  created. 

6 


Zjsia,  m)  -  0.5  X 


1  + 


i=l 


(48) 
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where 


(p\{d)  =  (1  -  of  +  6fl(l  -  of , 

<P2(a)  -  15fl^(l  - 

<P3(a)  -  20fl^(l  -  a)^, 

<P4(a)  =  -  a)^ 

(psia)  -  6fl^(l  -  a), 

<P6(a)  -  fl®, 

where  a  =  x/L  is  the  normalized  x-coordinate,/3,((y)  are  Bezier  coefficients,  which  are 
i.i.d.  randomly  sampled  from  the  uniform  distribution  'Z/(-0.1, 0.1).  At  the  beginning, 
we  assume  all  the  microstructures  to  be  identical  (same  grain  sizes  and  texture).  The 
only  difference  between  ingot  samples  is  the  random  shape  of  the  upper  surface.  All 
raw  ingots  are  then  used  as  an  input  to  a  deterministic  flat-die  forging  process,  during 
which,  their  wavy  surfaces  are  flattened  under  strain  rate  v  =  0.01  (Fig.  4).  Since 
all  workpieces  go  through  distinct  deformation  processes  due  to  their  unique  surface 
shapes,  the  resultant  microstructures  will  vary  from  point  to  point  and  from  sample  to 
sample.  The  resultant  microstructures  after  pre-processing  are  collected  as  the  data¬ 
base  of  the  following  stochastic  simulation.  They  will  be  adopted  to  build  the  reduced- 
order  model.  In  the  next  step  stochastic  simulation,  new  microstructure  samples  will 
be  sampled  and  assigned  to  a  regular-shaped  workpiece.  The  flattened  workpieces  in 
the  pre-process  are  abandoned.  The  information  about  the  generation  of  the  random 
microstructure  samples  is  totally  blind  to  the  following  stochastic  forging  simulation. 

Since  the  crystal  plasticity  constitutive  model  adopted  here  only  updates  grain  ori¬ 
entations  while  leaving  grain  size  untouched,  the  uncertainty  source  of  the  stochastic 
simulation  is  the  texture  of  microstructures,  which  has  been  proven  to  have  great  effect 
on  the  mechanical  response  and  properties  of  polycrystals  [12].  The  dimensionality 
of  the  input  is  x  x  n feature  —  60  X  4  X  60  =  14400,  where  is  the  num¬ 
ber  of  elements  in  the  macroscale  discretization,  is  the  number  of  Gauss  points  of 
an  element,  and  n feature  is  the  dimensionality  of  the  random  feature  that  describes  a 
microstructures  (20  set  of  3-dimensional  Rodrigues  parameters  in  here).  In  what  fol- 
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Figure  4:  Left:  Initial  ingot  with  random  upper  surface  and  identical  microstructures.  Right:  Flattened  ingot 
having  various  resultant  microstructures. 


lows,  the  vector  A(x,  s,  at)  is  written  as  a  scalar  A(x,  s,  a>)  since  the  only  random  source 
is  grain  orientation.  The  number  of  features  that  describe  a  microstructure  is  there¬ 
fore  being  m  -  I,  leading  to  the  fact  that  A(x,  s,  (o)  is  a  scalar  for  a  specific  x,  s  and 
CO.  Consequently,  the  reduced  macroscale  modes  {<t),(x,  and  reduced  represen¬ 
tations  /  =  I,. .  .,d  are  also  scalars  for  a  specific  x  and  co.  We  will  adopt 

the  aforementioned  two-step  KLE  to  reduce  the  dimensionality  of  the  stochastic  input 
space  driven  by  the  1000  sets  of  microstructure  samples.  The  information  of  how  these 
samples  are  generated  is  blind  to  the  model  reduction  process.  The  reduced  random 
variables  will  be  mapped  to  standard  multivariate  uniform  distribution  ('Z/(-l,  1))  fol¬ 
lowing  the  PC  expansion  through  the  non-intrusive  projection.  New  samples  will  be 
drawn  from  the  reduced  space  and  reconstructed  to  be  the  input  to  the  multiscale  phys¬ 
ical  simulator.  Monte  Carlo  simulation  is  employed  to  solve  the  underlying  stochastic 
equations  in  conjunction  with  the  multiscale  deterministic  forging  solver. 

The  procedure  of  the  complete  uncertainty  quantification  is  illustrated  in  Fig.  5 
and  summarized  below,  (a)  to  (b);  Given  a  number  of  initial  ingot  samples,  compute 
the  separated  mesoscale  and  macro-random  coupled  modes  using  the  bi-orthogonal 
KLE.  (b)  to  (c):  Project  the  macro-random  coupled  modes  to  low-dimensional  space 
through  a  second-level  KLE.  (c)  to  (d)  Map  the  reduced  stochastic  space  to  a  known 
(e.g.  uniform)  distribution  using  PCE.  (d)  to  (c):  Generate  new  samples  in  the  known 
low-dimensional  distribution,  and  find  their  counterparts  in  the  reduced  surrogate  space 
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Bi-orthogonal  KLE 


I 

_  d 

A(x,  s,  w)=5  A(x,  s)  +  ^  ^/A 'i',  (s)® ,  (*.  ®) 


(a)  Multiscale  microstructure  sample 

Forging 

v.,-(P,.)*+f,.=o 

T-r„  =  aMb^ 


Stochastic  Multiscale 
Simulation 


(b)  Macro-rancfom  coupled  modes 
~|  Second-level  KLE 

<l».(x,0))*O  (x) 


(e)  Properties  distribution 


(d)  Uniform  distribution 


(c)  Reduced  surrogate  space 


Figure  5:  Procedure  of  the  stochastic  multiscale  simulation  for  quantifying  variability  of  properties  of  forging 
disks  due  to  microstructure  uncertainty. 


through  PCE.  (c)  to  (b)  Recover  macroscale  modes  via  KLE.  (b)  to  (a):  Reconstruct 
physical  representation  (microstructures  of  the  workpiece)  of  new  samples  using  bi- 
orthogonal  KLE  after  obtaining  the  macroscale  modes,  (a)  to  (e):  Perform  multiscale 
forging  simulations  to  obtain  the  properties  of  the  reconstructed  samples.  Repeating 
(d)-(c)-(b)-(a)-(e)  multiple  times,  the  statistics  of  properties  of  the  products  can  be 
evaluated. 


6.  Numerical  Examples 

We  will  next  validate  the  bi-orthogonal  decomposition  strategy  for  reducing  the 
complexity  of  stochastic  multiscale  input.  Examples  comparing  reconstructed  features 
and  initial  sample  features,  as  well  as  their  corresponding  properties  after  forging,  are 
demonstrated. 
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6.1.  Construction  and  validation  of  the  reduced-order  model 

As  described  in  the  previous  section,  1000  sets  of  preforms  with  correlated  mi¬ 
crostructures  (textures)  that  resulted  from  the  same  preprocess  are  the  input  to  the 
stochastic  simulation.  The  variability  of  mechanical  properties  of  the  products  is  in¬ 
duced  by  the  randomness  of  these  initial  microstructures.  Since  each  macropoint  of 
the  workpiece  associates  with  a  20-grain  microstructure,  the  total  dimensionality  of  the 
input  is  14400  according  to  the  calculation  in  the  previous  section.  It  is  impossible  to 
explore  such  a  high  dimensional  space  and  investigate  the  variability  of  corresponding 
properties  of  products.  The  bi-orthogonal  decomposition,  followed  by  a  second-level 
KLE,  is  applied  to  the  1000  14400-dimensional  samples.  We  will  determine  the  opti¬ 
mal  dimensionality  of  the  reduced  space  by  the  proportion  of  energy  captured  by  the 
first  few  principal  components. 

In  bi-orthogonal  KLE,  the  random  energy  of  the  k-th  macro-random  coupled  mode 
is  defined  by  [19] 

Ef  \cj)  :  =  pf ’of >(x,  tu)  ■  of >(x,  co)dx,  (49) 

The  expectation  of  the  random  energy  is  therefore 

£f’  -  (£f’>-  (50) 


Since  the  mean  feature  A  defined  by  Eq.  (28)  has  been  extracted  from  the  centerized 
samples  A,  which  are  the  data  base  of  the  expansion,  the  inner  product  Eq.  (30)  defined 
with  h  -  Q  should  be  0.  As  a  result,  h  -  0  is  not  appropriate  for  constructing  the 
covariance  matrix  and  performing  model  reduction.  The  covariance  matrix  when  h  =  1 
is  identical  with  the  one  when  h  —  2,  due  to  Eq.  (37).  Therefore,  we  select  h  =  1  as  the 
inner  product  in  the  macroscale  for  the  current  work. 

Computing  the  energy  expectation  for  all  spatial  modes,  the  energy  spectrum  is 
plotted  in  Eig.  6  for  h  -  \.  The  energy  proportion  captured  by  the  first  d  modes  is 
defined  as 


jQi) 

Energy 


(d)  = 


X'M  pW  ’ 


where  the  mean  energies  are  sorted  in  descending  order. 


(51) 
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Figure  6:  The  energy  and  eigenvalue  spectmm  of  initial  samples.  The  value  of  y-axis  is  the  total  energy 
proportion  captured  by  the  first  x  principal  components. 


It  is  observed  that  the  first  few  energy  components  capture  most  of  the  total  energy. 
The  energy  spectrum  overlaps  with  the  eigenvalue  spectrum  defined  by 


’jih) 
Eigenvalue 


(d)^ 


Y'M  (h) 

^J=\Pi 


(52) 


To  effectively  reduce  the  complexity  while  preserving  most  of  the  features  of  the 
initial  samples,  we  truncate  the  bi-orthogonal  KLE  expansion  keeping  only  the  first 
3  modes,  which  captures  almost  95%  of  the  total  energy.  This  truncation  strategy  is 
similar  to  that  in  conventional  KLE. 

Remark  1:  Note  that  if  the  mean  feature  is  defined  as  [19] 


f  <A(x,  s,  6>))t/s,  (53) 

\M\  Jm 

the  bi-orthogonal  decomposition  gives  different  results.  In  this  case  the  first  type  of 
macroscale  inner  product  is  not  zero  any  more,  and  the  second  type  covariance  €*^^^(8,  s) 
becomes  small.  The  micro  and  macro  modes  become  different  from  the  ones  presented 
earlier. 

Remark  2:  The  overlap  of  the  eigenvalue  and  energy  spectrums  in  Eig.  6  is  a  result  of 
the  current  problem  setup.  Indeed,  combining  Eqs.  (49)  and  (50),  we  obtain 

=  r  Pf  «’(X’ ")  ■  (X- "))  dx.  (54) 

When  h  —  I,  and  considering  the  orthogonality,  we  obtain 
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The  energy  level  of  the  k-th  mode,  in  general,  does  not  reflect  the  magnitude  of  the 
corresponding  eigenvalue.  In  another  words,  a  mode  that  possesses  large  energy  does 
not  necessarily  correspond  to  a  large  eigenvalue.  An  example  is  shown  in  Fig.  7  for 
the  case  that  A  is  defined  by  Eq.  (53)  and  h  -  2.  This  outcome  is  consistent  with  the 
results  provided  in  [19].  These  facts  complicate  the  performance  of  the  bi-orthogonal 
decomposition.  Therefore,  alternative  setups  [Xah  -  1  and  Eq.  (28))  are  not  used  in  the 
current  work. 


Figure  7:  The  eigenvalue  and  mean  energy  of  the  bi-orthogonal  KLE.  Both  quantities  are  normalized  with 
respect  to  their  maximum  values. 

The  reduced  spatial  modes  jOi,  <1)2,  <1>3)  are  presented  by  1000  realizations  of  240- 
dimensional  vectors,  because  the  mesoscale  modes  have  been  segregated.  The  next 
task  is  to  separate  the  random  variables  from  spatial  coordinates  using  a  second-level 
KLE,  which  results  in  further  reduction  of  the  random  space.  Eor  each  <&,  {i  =1,2, 3), 
we  perform  an  independent  KLE  and  keep  the  largest  r,  components  that  capture  more 
than  95%  of  the  total  energy  of  <!>,■.  The  energy  spectrum  of  the  three  modes  are  plotted 
in  Eig.  8.  The  number  of  preserved  components  are  ri  -2,r2-  3,  r^,  -  3,  respectively. 
The  dimensionality  of  the  final  reduced  space  is  therefore  r  -  ri  +  r2  +  r^,  =8. 
Remark  3:  It  is  interesting  to  note  that  the  differences  between  large  eigenvalues  and 
small  eigenvalues  of  the  macro-modes  <&,  reduces  as  i  increases.  To  capture  95%  of  the 
total  energy,  only  the  largest  2  eigenvalues  is  enough  for  <I)i,  while  for  ^2  and  <1)3  three 
eigenvalues  are  needed.  We  also  examined  the  macro-modes  that  correspond  to  lower 
energy  in  the  bi-orthogonal  KLE  and  discovered  that  7,  9,  35,  and  108  principle  com¬ 
ponents  are  needed  to  capture  95%  energy  of  macro  modes  from  <1)4  to  (by,  respectively. 
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Figure  8:  The  eigenvalue  spectrum  of  three  macro-random  coupled  modes.  The  value  of  }j-axis  is  the  total 
energy  proportion  captured  by  the  first  x  principal  components.  Only  the  first  50  dimensions  are  shown. 

The  increase  of  the  dimensions  (e.g.  number  of  principal  components)  that  are  neces¬ 
sary  to  capture  the  same  proportion  of  the  total  energy  is  dramatic,  when  the  energy 
captured  by  the  macro-modes  O,-  decreases.  For  this  reason,  keeping  a  small  number  of 
<t);  is  of  great  importance  to  reduce  the  stochastic  space.  In  the  curi'ent  example,  d  -  3 
is  the  optimal  choice. 

As  we  discussed  before,  the  reduced  surrogate  space  of  input  after  two-step  KLE 
needs  to  be  projected  to  a  well-shaped  equivalent  space  through  PCE.  Therefore,  new 
samples  can  be  easily  drawn.  Legendre  polynomials  are  selected  as  we  will  map  the 
reduced-order  random  variables  to  uniform  distributions  'Hi- 1,1).  The  order  of  the  PC 
basis  is  set  to  be  12,  which  gives  accurate  estimation  to  the  distributions  of  the  reduced 
representations.  We  plot  and  compare  the  PDEs  of  initial  reduced  representations  and 
new  samples  in  Pigs.  9-11.  The  distributions  of  the  initial  reduced  representations 
(pi,i  —  .  ,d,  i  -  .  ,ri  are  computed  from  the  histogram  of  the  given  1000  initial 

samples  derived  by  the  two-step  KLE.  On  the  other  hand,  10000  new  samples  are 
randomly  sampled  from  the  uniform  distribution  and  mapped  to  the  suri'ogate  space 
via  PCE.  A  great  consistence  of  the  two  curves  is  observed  and  thus  the  PCE  on  reduced 
random  variables  is  successful. 

To  check  the  performance  of  the  multiscale  model  reduction,  we  compare  a  recon¬ 
structed  sample  with  its  original  realization.  An  14400-dimensional  array.  A,  contain¬ 
ing  the  texture  information  of  all  initial  microstructures  over  a  workpiece  realization 
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Figure  9:  Mai'ginal  PDFs  of  initial  low-dimensional  representations  (i.e.  random  variables),  =  1,2, 


corresponding  to  the  first  spatial  mode  <I>i  (the  reduced  representations  obtained  after  two-step  KLE  on  the 
1000  given  texture  samples)  and  identified  random  variables  obtained  using  PCE  (reconstructed  through  PCE 
on  10000  randomly  generated  samples  from  uniform  distribution).  The  distributions  are  constructed  through 
kernel  density  based  on  data. 

is  projected  to  the  8-dimensional  uniform  distribution  through  the  two-step  KLE-PCE 
process.  These  8  reduced  variables  within  (-1, 1)  are  then  mapped  back  to  a  14400- 
dimensional  texture  array  through  the  inverse  PCE  and  KLE.  We  hrst  compare  the  re¬ 
stored  spatial  modes  from  the  reduced  variables  with  the  initial  modes  obtained  through 
bi-orthogonal  KLE  on  the  original  sample.  The  3  modes  capturing  most  energy  are 
shown.  It  is  observed  in  Eig.  12  that  the  restored  are  close  to  the  original  ones. 

We  further  reconstruct  the  texture  realization  in  the  physical  space  based  on  the 
restored  spatial  modes  through  Eq.  (38).  The  restored  texture  is  compared  with  the 
original  one  and  a  relative  error  dehned  as 


M  Ai  _  Ai 


liestored 


(55) 


is  computed.  In  Eq.  (55),  M  is  the  dimensionality  of  the  texture  array  of  the  entire 
workpiece  (here,  M  —  14400),  and  are  the  Rodrigues  parameters  of 

the  original  and  restored  samples,  respectively.  Notice  that  each  orientation  is  described 
by  3  Rodrigues  parameters.  In  the  current  setup,  we  put  them  all  in  an  1 -dimensional 
array.  Eigure  13  (a)  shows  the  reconstructed  and  original  texture  of  the  microstructure 
associated  with  one  macropoint  (for  a  single  microstructure).  The  relationship  between 
the  two  samples  throughout  the  entire  workpiece  is  depicted  in  Eig.  13(b).  We  observe 
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Figure  10:  Marginal  PDFs  of  initial  low-dimensional  representations  (i.e.  random  variables),  i  =  1. 2, 3, 
corresponding  to  the  second  spatial  mode  <D2  (the  reduced  representations  obtained  after  two-step  KLE  on 
the  1000  given  texture  samples)  and  identified  random  variables  obtained  using  PCE  (reconstructed  through 
PCE  on  10000  randomly  generated  samples  from  uniform  distribution).  The  distributions  are  constructed 
through  kernel  density  based  on  data. 

a  nearly  straight  line  in  the  original-restored  texture  plot,  meaning  that  the  two  samples 
are  almost  identical.  The  deviation  of  the  restored  samples  from  the  original  one  is 
small.  The  great  consistence  of  the  reconstructed  and  original  samples  is  obtained. 
The  relative  error  is  ~  4.26%. 

6.2.  Stochastic  multiscale  forging  simulation 

After  establishing  the  connection  between  the  microstructure  space  and  the  reduced 
suiTogate  space,  we  are  ready  to  draw  random  samples  for  the  investigation  of  the 
mechanical  properties  of  workpieces  whose  microstructures  are  statistically  similar  to 
the  given  data.  The  mean  and  standard  deviation  of  the  equivalent  strain,  stress,  and 
strength  fields  of  the  forged  workpiece  are  of  interest.  The  effective  strain  and  stress  are 
evaluated  using  Eqs.  (21)  and  (22).  The  effective  strength  is  measured  by  the  average  of 
the  slip  resistances  over  all  slip  systems  of  all  grains  in  a  microstructure  (Eq.  (24)).  The 
mean  and  standard  deviation  fields  computed  based  on  4032  MC  samples  randomly 
generated  from  the  reduced-space  are  plotted  in  Eigs.  14  and  15.  The  fields  computed 
directly  from  the  1000  initial  samples  are  also  plotted  (in  the  same  figures)  in  com¬ 
parison  with  the  reconstructed  results.  The  mean  fields  of  properties  computed  from 
reconstructed  samples  are  close  to  the  ones  computed  from  the  initial  samples.  This 
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Figure  11:  Marginal  PDFs  of  initial  low-dimensional  representations  (i.e.  random  variables),  (pyj=  1, 2, 3, 
corresponding  to  the  third  spatial  mode  <I>3  (the  reduced  representations  obtained  after  two-step  KLE  on  the 
1000  given  texture  samples)  and  identified  random  variables  obtained  using  PCE  (reconstructed  through  PCE 
on  10000  randomly  generated  samples  from  uniform  distribution).  The  distributions  are  constructed  through 
kernel  density  based  on  data. 


is  consistent  with  the  bi-orthogonal  decomposition  setup.  The  standard  deviation  of 
the  properties  of  reconstructed  samples,  however,  shows  deviation  from  that  computed 
using  the  initial  samples.  This  is  because  the  limited  number  of  given  samples  are  not 
enough  to  represent  the  entire  random  microstructure  space  (especially  the  higher  order 
statistics).  On  the  other  hand,  the  reconstructed  samples  are  generated  from  the  surro¬ 
gate  space  which  is  built  to  efficiently  represent  the  complete  microstructure  space. 
Random  samples  from  the  reduced-order  model  reveal  features  that  cannot  be  captured 
by  the  given  initial  samples.  The  gained  efficiency  in  sampling  in  the  low-dimensional 
surrogate  microstructure  space  is  prominent. 

A  convergence  test  is  also  conducted  using  8064  random  MC  samples.  The  com¬ 
parison  of  the  mean  and  standard  deviation  between  4032  and  8064  samples  are  plotted 
in  Figs.  16  and  17.  The  relative  difference  of  quantities  between  the  two  sets  of  sim¬ 
ulations  defined  as  (Pgoef  -  P4032)  I Pm64,  where  Pff  is  the  quantity  evaluated  using  N 
MC  samples,  is  plotted  in  Fig.  18.  From  the  difference  we  see  that  the  mean  fields  of 
the  two  simulations  are  almost  the  same.  The  relative  error  of  standard  deviation  fields 
is  larger  than  that  of  the  mean  field.  The  largest  error  is  around  0.05. 

In  order  to  test  the  convergence  of  the  bi-orthogonal  decomposition  model  reduc¬ 
tion  scheme,  we  next  keep  more  components  in  the  second-level  KLE  so  that  they 
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Figure  12:  Comparison  of  reconstructed  and  initial  spatial  modes  of  a  single  texture  sample.  The  original 
modes  are  obtained  by  projecting  a  initial  texture  sample  to  eigenbasis  through  bi-orthogonal  KLE.  The 
reconstructed  modes  are  recovered  from  the  low-dimensional  representations  via  PCE  and  second-level  KLE. 
The  dimensionality  of  the  reduced  representations  of  Oi,  02,  and  O3  are  2,  3,  and  3,  respectively. 


Figure  13:  (a)  The  reconstructed  and  original  texture  of  the  microstructure  at  a  single  macropoint,  (b) 
The  comparison  between  the  reconstructed  texture  and  its  original  realization  of  the  entire  workpiece.  The 
reconstructed  texture  is  obtained  from  a  8-dimensional  representation. 


capture  99%  energy  of  the  macro-modes.  The  dimensionality  of  the  reduced  space  be¬ 
comes  r  -  ri  +  r2  +  r-i  =  3h-7-i-8  =  18.  The  reconstructed  macro-modes  are  obviously 
closer  to  their  original  samples  as  shown  in  Fig.  19.  Similarly,  the  reconstructed  texture 
has  smaller  error,  s  -  0.0398,  comparing  with  the  original  sample  (Fig.  20). 

The  mean  and  standard  deviation  of  effective  strain,  stress  and  strength  fields  are 
plotted  in  Fig.  21.  The  relative  difference  of  fields  defined  as  (Pig  -  Pg)/Pi8,  where 
Pd  is  the  quantity  evaluated  from  cZ-dimensional  reduced  space,  is  shown  in  Fig.  22. 
It  is  observed  that  keeping  18  reduced  variables  gives  very  similar  results  as  keeping 
8  low-dimensional  representations,  since  the  total  energy  captured  by  the  two  cases  is 
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Figure  14:  Mean  field  of  the  properties  of  the  forged  product.  Upper:  results  extracted  from  1000  initial 
samples;  lower:  results  evaluated  through  4032  MC  samples  randomly  generated  from  the  8-dimensional 
reduced  space,  (a)  elfective  strain,  (b)  elfective  stress,  (c)  effective  strength. 


close.  The  number  of  samples  used  here  is  8064. 

The  distributions  of  properties  of  any  point  on  the  solid  can  also  be  computed.  In 
Fig.  23,  we  plot  the  equivalent  strain,  stress,  and  strength  distributions,  as  well  as  the 
convex  hull  of  these  three  quantities,  at  a  single  spot  of  the  workpiece,  where  the  equiv¬ 
alent  strain  is  large.  All  distributions  and  the  convex  hull  [34]  are  evaluated  according 
to  the  results  of  4032  randomly  generated  samples  from  the  8-dimensional  reduced 
space  in  the  MC  simulation  just  discussed. 


7.  Conclusions 

A  multiscale  model  reduction  scheme  based  on  the  bi-orthogonal  KLE,  conven¬ 
tional  KLE  and  PCE  was  presented.  The  basic  idea  is  to  decompose  the  multiscale 
random  field  into  a  few  orthogonal  modes  in  different  (macro  and  meso)  scales  and 
separate  the  inherent  random  variable  from  the  two  scale  coordinates.  A  non-intrusive 
projection  strategy  is  employed  to  map  the  reduced  representations  after  the  two-step 


31 


Initial  samples 


(a)  (b)  (c) 


Figure  15:  Standai'd  deviation  field  of  the  properties  of  the  forged  product.  Upper:  results  extracted  from 
1000  initial  samples;  lower:  results  evaluated  through  4032  MC  samples  randomly  generated  from  the  8- 
dimensional  reduced  space,  (a)  effective  strain,  (b)  effective  stress,  (c)  effective  strength. 

KLE  to  a  multivariate  uniform  distribution.  The  reconstructed  realizations  show  agree¬ 
ment  with  the  initial  microstructure  samples  that  are  given  as  the  known  information. 

In  the  context  of  polycrystalline  processes,  the  multiscale  random  field  is  location- 
dependent  high-dimensional  random  microstructure  features,  which  is  reduced  to  a 
low-dimensional  surrogate  space.  By  sampling  in  the  reduced  surrogate  space,  we  can 
equivalently  exploit  the  original  high-dimensional  microstructure  space.  Properties  of 
a  continuum  workpiece  subjected  to  forging  are  evaluated  by  a  multiscale  solver  which 
couples  FE  large  deformation  simulator  with  crystal  plasticity  constitutive  model.  The 
mean  and  standard  deviation  of  the  equivalent  strain,  stress,  and  strength  of  the  fi¬ 
nal  product  are  computed  using  MC  and  ASGC  methods.  It  is  seen  that  the  reduced 
model  captures  most  features  of  the  full  model  making  feasible  to  perform  large  scale 
stochastic  multiscale  simulation.  Future  studies  will  focus  on  the  model  reduction  of 
realistic  microstructures  described  by  pixels  rather  than  statistical  features  as  in  the 
present  work.  Moreover,  a  more  robust  strategy  of  decomposing  macro-random  modes 
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4032  samples 


(a)  (b)  (c) 


Figure  16:  Convergence  test  of  the  mean  field  of  the  properties  of  the  forged  product.  Upper:  results  extracted 
from  4032  MC  samples  randomly  generated  from  the  8-dimensional  reduced  space;  lower:  results  evaluated 
through  8064  MC  samples  randomly  generated  from  the  8-dimensional  reduced  space,  (a)  effective  strain, 
(b)  effective  stress,  (c)  effective  strength. 

after  bi-orthogonal  KLE  is  of  interest. 
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Figure  17:  Convergence  test  of  the  standard  deviation  field  of  the  properties  of  the  forged  product.  Upper: 
results  extracted  from  4032  MC  samples  randomly  generated  from  the  8-dimensional  reduced  space;  lower: 
results  evaluated  through  8064  MC  samples  randomly  generated  from  the  8-dimensional  reduced  space,  (a) 
effective  strain,  (b)  effective  stress,  (c)  effective  strength. 
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Figure  18:  Relative  difference  of  the  mean  and  standard  deviation  field  of  the  properties  of  the  forged  product 
computed  by  8064  and  4032  MC  samples  drawn  from  the  8-dimensional  reduced  space.  Upper:  difference 
of  mean  fields;  lower:  difference  of  standard  deviation  fields,  (a)  effective  strain,  (b)  effective  stress,  (c) 
effective  strength. 


Figure  19:  Comparison  of  reconstructed  and  initial  spatial  modes  of  a  single  texture  sample.  The  original 
modes  are  obtained  by  projecting  a  initial  texture  sample  to  eigenbasis  through  bi-orthogonal  KLE.  The 
reconstructed  modes  are  recovered  from  the  low-dimensional  representations,  with  higher  dimensions,  via 
PCE  and  second-level  KLE.  The  dimensionality  of  the  reduced  representations  of  Oi,  <1)2,  and  <1)3  are  3,  7, 
and  8,  respectively. 
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Figure  20:  (a)  The  reconstructed  and  original  texture  of  the  microstnacture  at  a  single  macropoint,  (b) 
The  comparison  between  the  reconstructed  texture  and  its  original  realization  of  the  entire  workpiece.  The 
reconstructed  texture  is  obtained  from  a  18-dimensional  representation. 


Figure  21:  The  mean  and  standard  deviation  fields  of  effective  strain,  stress,  and  strength  computed  based  on 
random  microstructures  reconstructed  from  18-dimensional  reduced-order  representations. 
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Figure  22:  The  relative  eiTor  of  mean  and  standard  deviation  fields  of  effective  strain,  stress,  and  strength 
computed  based  on  random  microstnactures  reconstructed  from  8-dimensional  and  18-dimensional  reduced- 
order  representations. 
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Figure  23:  Variability  of  properties  at  one  single  point  of  the  forged  disk  with  random  microstructures,  (a) 
A  convex  hull  showing  the  envelope  of  the  three  properties,  (b)  Equivalent  stress  distribution,  (c)  Equivalent 
strength  distribution,  (d)  Equivalent  strain  distribution. 
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